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ABSTRACT 


This  technical  report  sunroarizes  a  new  chemical  equilibrium  computer  code 
developed  for  general  reacting  systems,  and  presents  equilibrium  results  for 
several  reacting  pairs  in  liquid  metal  fuel  combustion.^  The  content  herein  is 
based  on  the  M.S.  dissertation  submitted  by  ,  Department  of 

MechanicaJ_£of4fleeriiig,  UniversTty^of  Wisconsin— Milwaukee. 


ASpecif ically,  the  object  of  this  report  is  to  develop  a  general  computer 
code  for  calculating  complex  chemical  equilibrium  of  nonideal,  multiphase, 
electrolytic  mixtures.  The  necessary  thermodynamic  foundation  is  provided  and 
thermodynamic  potentials  which  describe  the  equilibrium  state  of  a  system  are 
developed  from  the  fundamental  relation  using  Legendre  transforms.  For  a 
system  at  a  given  temperature  and  pressure,  the  Gibbs'  potential  function 
reaches  a  minimum  value  at  equilibrium.  The  minimization  of  this  function, 
subject  to  element  and  charge  conservation  constraints,  is  viewed  as  a 
constrained  optimization  problem,  which  is  solved  by  the  method  of  Lagrangian 
multipliers.  A  new  general  computer  code,  called  CEC-NMS,  has  been  developed 
and  is  presented  here  for  computing  chemical  equilibrta.  The  evaluation  and 
calculation  of  necessary  program  inputs  such  as  thermodynamic  data  and 
specifying  the  relative  amounts  of  reactants  are  discussed.  Equilibrium 
results,  including  temperature,  density,  and  species  concentrations  as  a 
function  of  mixture  fraction  are  presented  for  systems  using  the  following 
reactants: 


»^^g)-^»3(aq)-»2°(L) 

2.  Cl^^gj-Na^L) 
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Several  different  modeling  techniques  are  used  to  accurately  estimate  the 
activity  coefficients  of  the  species  in  the  above  systems.  The  liquid  phase 
of  the  electrolytic  solution  is  calculated  using  a  modified  Pitzer  formulation 
for  multicomponent  electrolytes  with  molecular  species.  The  gas  phase  of  the 
electrolytic  solution  is  modeled  using  a  pressure-explicit  second  order  virial 
equation.  The  pure-component  and  cross-component  second  virial  coefficients 
are  predicted  using  O'Connell's  generalized  method.  For  the  systems 
containing  lithium  and  sodium,  the  gas  phase  is  assumed  to  be  ideal,  and  the 
activity  coefficients  of  the  liquid  metal  and  molten  salt  comprising  the 
immiscible  liquid  phases  are  calculated  using  the  van  Laar  model. 
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1 .  INTRODUCTION 


1 .1  Definition  of  Equilibrium 

As  an  introduction  to  this  work,  it  is  useful  to  present  several 
fundamental  definitions.  Callen  (1985)  defines  thermodynamic  equilibrium  as 
follows: 

...  in  all  systems  there  is  a  tendency  to  evolve  toward  states  in 
which  the  properties  are  determined  by  Intrinsic  factors  and  not  by 
previously  applied  external  influences.  Such  simple  terminal  states 
are,  by  definition,  time  independent.  They  are  called  equilibrium 
states.  (Callen,  1985,  p.  13) 

In  addition,  he  emphasizes  the  importance  of  the  concept  of  equilibrium  by 
calling  the  determination  of  equilibrium  states  "...  the  single,  all- 
encompassing  problem  of  thermodynamics"  (Callen,  1985,  p.  26). 

The  thermodynamic  equilibrium  of  a  simple  system  (i.e.,  systems  that  are 
not  acted  on  by  electric,  magnetic,  or  gravitational  fields,  and  are 
macroscopically  homogeneous)  requires  three  distinct  kinds  of  equilibrium 
(Kuo.  1986,  p.  9): 

1.  Thermal  equilibrium  -  exists  when  all  parts  of  the  system  are  at  the 
same  temperature. 

2.  Mechanical  equilibrium  -  exists  when  there  are  no  unbalanced  forces 
within  the  system. 

3.  Chemical  equilibrium  -  exists  when  a  system  has  no  tendency  to 
undergo  a  spontaneous  change  In  chemical  composition. 

1 .2  Applications 

The  topic  of  this  report  is  the  computation  of  chemical  equilibrium  for 
multiphase  systems.  Numerous  applications  exist  for  chemical  equilibrium 
calculations.  Van  Zeggeren  and  Storey  (1970)  describe  four  major  areas: 


1.  Determining  the  specific  impulse  of  a  propellant  in  a  rocket  motor. 

2.  Calculation  of  the  properties  of  explosives. 

3.  Chemical  processing. 

4.  Calculating  the  behavior  of  multiphase  biological  cell  systems. 

In  addition,  Smith  and  Nissen  (1982)  discuss  applications  of  chemical 
equilibrium  analysis  in  chemical  kinetics,  inorganic  and  organic  chemistry, 
energy  conversion,  analytical  chemistry,  and  environmental  chemistry. 

The  impetus  for  this  particular  foray  into  the  field  of  chemical 
equilibrium  involves  modeling  the  reaction  zone  structure  of  a  turbulent, 
reacting  jet  discharging  into  a  stagnant  bath.  As  developed  by  Shearer  et  al. 
(1979),  Chen  and  Faeth  (1982,  1983)  and  Chan  et  al.  (1987,  1988),  the 
analytical  model  assumes  local  thermodynamic  equilibrium  at  each  point  in  the 
flow.  As  an  input  it  requires  the  calculation  of  equilibrium  compositions  for 
all  possible  ratios  of  oxidant  to  fuel.  This  model  allows  the  prediction  of 
the  concentration  of  chemical  species,  velocity,  temperature,  and  mixture 
density  at  all  points  in  the  flow  field.  It  has  been  used  with  reasonable 
success  by  this  laboratory  at  the  University  of  Wisconsin— Milwaukee  (UWM)  in 
the  past  (Chan  et  al.,  1987,  1988). 

1 .3  Computational  Methods 

Numerous  methods  have  been  proposed  for  computing  chemical  equilibrium  in 
the  last  40  years.  These  are  reviewed  by  Gautam  and  Seider  (1979),  Smith  and 
Missen  (1982),  van  Zeggeren  and  Storey  (1970),  and  Zeleznik  and  Gordon 
(1968).  While  there  are  several  ways  of  grouping  the  approaches,  Zeleznik  and 
Gordon  (1968)  divide  the  general  purpose  methods  into  the  equilibrium  constant 
approach  and  the  free  energy  minimization  approach.  In  the  former  method 
chemical  reaction  equations  are  written  to  describe  the  relevant  equilibria. 
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Equilibrium  constants  for  these  equations  are  calculated  which  relate  the 
concentrations  of  the  species  in  the  reaction.  The  resulting  system  of 
nonlinear  equations  may  be  solved  using  a  variety  of  techniques.  The  latter 
method  involves  a  search  for  the  minimum  in  the  Sibbs  free  energy  of  the 
system.  While  the  equilibrium  constant  method  as  formulated  by  Brinkley 
(1947)  and  the  free  energy  minimization  method  as  formulated  by  White  et  a1. 
(1958)  can  be  shown  to  be  computationally  equivalent  (Zeleznik  and  Gordon, 
1960;  Smith  and  Missen,  1982),  the  latter  method  eliminates  the  necessity  of 
writing  chemical  equations. 

1 .4  Existing  Computer  Programs 

Several  computer  programs  are  available  for  the  solution  of  general 
chemical  equilibria.  Probably  the  best  known  is  CEC-72  (Gordon  and  McBride, 
1971).  This  computer  code  can  be  used  for  ideal,  multiphase  calculations, 
provided  the  total  mass  of  the  condensed  (i.e.,  liquid  and  solid)  phases  are 
no  more  than  a  few  percent  of  the  total  mass  of  the  system.  No  provision  is 
made  for  the  use  of  activity  coefficients  to  describe  nonideal  mixtures.  The 
thermodynamic  state  of  the  system  may  be  specified  by  assigning  temperature 
and  pressure  (T,P),  enthalpy  and  Pressure  (H,P),  entropy  and  pressure  (S,P), 
temperature  and  volume  (T,V),  internal  energy  and  volume  (U,V),  or  entropy  and 
volume  (S,V).  A  data  file  containing  thermodynamic  property  data  for  more 
than  400  chemical  species  over  a  temperature  range  of  300  to  5000  K  is 
provided. 

Another  popular  computer  program  is  SOLGASMIX  (Ericksson  1971,  1973,  and 
1975).  This  program  can  be  used  to  calculate  equilibria  for  nonideal, 
multiphase  systems,  provided  that  the  activities  of  the  species  can  be 
specified  on  a  mole  fraction  scale  (i.e.,  -»  1  as  1 ,  where 
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is  the  activity  coefficient  and  is  the  mole  fraction  of  species  i 
in  the  phase).  Therefore,  it  is  not  applicable  to  electrolytic  solutions, 
where  the  activities  (Jf  the  solutes  are  normally  based  on  molality.  The 
thermodynamic  state  may  be  specified  by  temperature  and  pressure  only.  All 
thermodynamic  data  must  be  supplied  by  the  user.  An  extension  of  this  program 
by  Besmann  (1977)  uses  the  ideal  gas  law  to  allow  the  calculation  of 
equilibria  at  constant  total  gas  volume.  Therefore,  the  thermodynamic  state 
may  also  be  specified  by  temperature  and  volume  (T,V). 

Ericksson  (1979)  also  developed  a  program  called  SOLGASWATER,  which  is 
designed  to  be  used  to  calculate  multiphase,  aqueous  equilibria.  The  system 
is  assumed  to  have  a  solvent  of  unit  activity  and  a  constant  volume  gas 
phase.  Relationships  specifying  the  activity  coefficients  of  the  solute 
species  must  be  provided  by  the  user.  Because  there  is  no  provision  for 
considering  nonideal  solvents,  accurate  predictions  of  equilibrium 
compositions  are  limited  to  dilute  solutions. 

Smith  and  Missen  (1982)  present  two  algorithms,  BNR  and  VCS,  useful  for 
calculating  ideal,  multiphase  equilibria  at  a  specified  temperature  and 
pressure.  The  BNR  program  is  based  on  the  Rand  method  (White  et  al.,  195B), 
and  the  VCS  program  is  based  on  a  formulation  using  extent  of  reaction 
variables.  Thermodynamic  data  must  be  supplied  by  the  user  for  both  programs. 

Zemaitis  et  al.  (19B6)  discuss  several  computer  programs  that  can  be  used 
for  electrolytic  solutions.  Generally,  the  publically  available  programs  are 
somewhat  limited  in  terms  of  features  and  applicability.  Several  proprietary 
programs,  however,  appear  to  be  quite  complete  and  powerful. 

1 .5  Objectives 

Work  in  this  laboratory  requires  a  computer  code  with  a  number  of 
features  not  found  (to  the  authors'  knowledge)  on  existing,  puolically 
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available  programs.  These  Include  the  ability  to  calculate  equilibria  of 
electrolytic  and  non-electrolytic,  multiphase,  nonideal  mixtures  at  either 
specified  temperature  and  pressure  (T,P)  or  specified  AH  and  pressure 
(AH,P)  (i.e.,  the  temperature  is  determined  indirectly  by  inclusion  of  an 
energy  equation  which  contains  a  specified  enthalpy  change  term  -  for  the 
adiabatic  case,  AH  =  0).  Flexibility  in  the  specification  of  activity 
coefficients  is  required,  so  that  they  may  be  written  in  whatever  convention 
(molality  or  mole  fraction)  is  relevant  to  the  particular  problem. 

This  report  presents  the  development  of  such  a  computer  code  which 
utilizes  a  variation  of  the  direct  Gibbs  free  energy  minimization  approach 
first  presented  by  White  et  al.  (1958)  while  working  at  Rand  Corporation. 
Section  2  attempts  to  provide  the  thermodynamic  framework  on  which  a 
discussion  of  chemical  equilibria  can  be  based.  Section  3  presents  the 
analytical  derivation  of  the  direct  minimization  approach.  Section  4 
describes  the  numerical  technique  used  to  solve  the  resulting  nonlinear 
equations.  The  computer  code  is  described  in  Section  5  and  several 
applications  of  the  program  are  demonstrated  in  Sections  6-9. 
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2.  THERMODYNAMIC  FUNDAMENTALS 


Callen  (1985)  bases  his  development  of  thermodynamic  theory  on  four 
postulates.  We  present  the  first  three  of  his  postulates  below.  The  fourth 
postulate  is  a  restatement  of  the  third  law  of  thermodynamics  and  is  not 
relevant  to  this  discussion. 


2.1  Postulates 

The  first  postulate  in  Callens'  development  proposes  the  existence  of  the 
equilibrium  state  (recall  the  definition  of  a  simple  system  given  in  Chapter 
1): 


Postulate  I.  "There  exist  particular  states,  (called  equilibrium 
states)  of  simple  systems  that,  macroscopically,  are  characterized 
completely  by  the  internal  energy  U,  the  volume  V,  and  the  mole 
numbers  N^ ,  M2,  ...,  M^  of  the  chemical  components."  (Callen, 

1985,  p.  13) 

Several  definitions  are  necessary  before  proceeding  with  the  next 
postulate.  We  define  a  composite  system  as  two  or  more  simple  systems.  We 
further  define  extensive  parameters  as  those  parameters  that  have  values  in  a 
composite  system  equal  to  the  sum  of  their  values  in  each  of  the  subsystems. 
The  second  postulate  proposes  the  existence  of  entropy  and  the  entropy  maximum 
principle: 

Postulate  II.  "There  exists  a  function  (called  the  entropy  S)  of 
the  extensive  parameters  of  any  composite  system,  defined  for  all 
equilibrium  states  and  having  the  following  property:  The  values 
assumed  by  the  extensive  parameters  in  the  absence  of  an  internal 
constraint  are  those  that  maximize  the  entropy  over  the  tnanifold  of 
constrained  equilibrium  states."  (Callen,  1985,  p.  27) 

The  relation  that  gives  the  entropy  as  a  function  of  the  extensive 

parameters  is  called  the  fundamental  relation.  If  the  fundamental  relation  of 

a  particular  system  is  known,  all  possible  thermodynamic  information  about  the 

system  can  be  derived  from  it.  The  fundamental  relation  in  entropic  form  can 

be  written  as: 
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S  =  S(U.  V.  N^,  . N^)  (2  1) 

The  third  postulate  specifies  certain  properties  of  entropy: 

Postulate  III.  "The  entropy  of  a  composite  system  is  additive  over 
the  constituent  subsystems.  The  entropy  is  continuous  and 
differentiable  and  is  a  monotonically  increasing  function  of  the 
energy."  (Callen,  1985,  p.  28) 

The  additivity  property  of  entropy  requires  that  the  entropy  of  a  simple 
system  be  a  homogenous  first-order  function  of  the  extensive  parameters.  That 
is,  if  all  the  extensive  parameters  of  a  system  are  multiplied  by  a  constant 
o,  the  entropy  is  multiplied  by  this  same  constant: 

S(aU.  aV,  aN^.  ...,  aN^)  =  aS(U,  V.  . N^^)  (2.2) 


2.2  Equations  of  State 

The  significance  of  the  differentiability,  continuity,  and  monotonic 

properties  of  the  entropy  are  that  they  allow  the  entropy  function  to  be 

inverted  with  respect  to  the  internal  energy.  The  internal  energy  is  then  a 

single-valued,  continuous,  and  differentiable  function  of  S,  V,  N^,  ..., 

M  .  The  fundamental  equation  can  then  be  written  in  the  form: 
n 

u  =  U(S,  V,  . N^)  (2.3) 

Computing  the  first  differential: 


dU  =  (3U/aS),.  „  dS 
V,N 


(3U/3V)^,N  dV 


n 

+  I 

i=l 


(2.4) 


The  partial  derivatives  in  the  preceding  equation  occur  so  frequently  that 
special  symbols  are  used  to  represent  them: 


(3U/3S)^  E  T 

(2.5) 

-(aU/3V)^  E  P 

(2.6) 
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=  “1 

(2.7) 

T,  P,  and  are  the  temperature,  pressure, 

1^^  component,  respectively. 

and  chemical  potential  of  the 

The  functional  relationships 

T  =  T(S,  V.  . N^) 

(2.8) 

P  =  P(S.  V,  . N^) 

(2.9) 

=  P^(S,  V.  . N^) 

(2.10) 

are  called  equations  of  state.  Since  the  fundamental  equation  Is  homogeneous 
first  order,  the  equations  of  state  are  homogeneous  zero  order.  This  means 
that  If  we  multiply  each  of  the  Independent  extensive  parameters  by  a  scalar 
a,  the  function  Is  left  unchanged: 

T(«S.  aV.  «N^ . aN^)  =  T(S,  V.  . N^)  (2.11) 

Parameters  for  which  this  is  true  are  said  to  be  intensive  parameters. 

2.3  Thermodynamic  Potentials 

A  mathematical  technique  known  as  the  Legendre  transformation  (Callen, 
1985,  pp.  137-151)  Is  performed  on  the  energy  representation  of  the 
fundamental  relation  (equation  2.3)  in  order  to  develop  fundamental  relations 
In  more  convenient  forms  for  particular  problems.  The  Legendre  transformed 
functions  are  called  thermodynamic  potentials. 

The  partial  Legendre  transform  of  U  that  replaces  the  entropy  by  the 
temperature  as  the  independent  variable  is  called  the  Helmholtz  potential  (F) 
or  Helmholtz  free  energy: 

F  =  F(T.  V.  N^,  . N^)  (2.12) 
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also 


F  =  U  -  TS  (2.13) 

As  a  result  of  the  entropy  maximum  principle  stated  in  postulate  II,  it  can  be 
shown  that  the  Helmholtz  potential  of  a  system  at  constant  temperature  is 
minimized  at  equilibrium. 

The  partial  Legendre  transform  of  U  that  replaces  the  volume  by  the 
pressure  as  an  independent  variable  is  called  the  enthalpy  (H): 

H  =  H(S.  P.  N^,  N2 . N^)  (2.14) 

also 

H  =  U  +  PV  (2.15) 

The  enthalpy  of  a  system  at  constant  pressure  is  minimized  at  equilibrium. 

The  partial  Legendre  transform  of  U  that  simultaneously  replaces  the 
entropy  by  the  temperature  and  the  volume  by  the  pressure  as  independent 
variables  is  called  the  Gibbs  potential  (G)  or  Gibbs  free  energy: 

G  =  G(T,  P,  N,,  . N^)  (2.16) 

also 

G  =  U  -  TS  +  PV  (2.17) 

The  equilibrium  state  of  a  system  at  constant  temperature  and  pressure  is  the 
state  in  which  the  Gibbs  potential  is  minimized. 

The  homogeneous  first-order  property  of  the  fundamental  relation 
described  above  permits  equation  2.3  to  be  written  in  a  different  form,  known 
as  the  Euler  relation  (Callen,  1985,  p.  59): 

U  =  TS  -  PV  +  +  ...  f  (2.18) 
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Substituting  the  Euler  relation  Into  equation  2.17,  we  get 


or,  In  more  compact  notation 


G  =  X  y^N. 
1=1  ’ 


(2.19) 


Therefore  the  Gibbs  potential  for  a  multicomponent  system  is  related  to  the 
chemical  potentials  of  the  individual  components. 
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3.  FORMULATION  OF  THE  DIRECT  MINIMIZATION  TECHNIQUE 


As  shown  in  Section  2,  the  condition  for  chemical  equilibrium  at  constant 
temperature  and  pressure  is  the  minimization  of  the  Gibbs  potential.  This 
minimization  is  subject  to  the  element  abundance  constraints  (conservation  of 
chemical  elements  making  up  the  species  of  the  system). 

3.1  Constraint  Equations 

The  constraining  element  abundance  equations  may  be  written  as 
n 

^i'^i  "  h’  k  =  1,  2,  ...,  m 

or,  equivalently 
n 

=  0;  k  *  1,  2 . m  (3.1) 

where  Aj^^  is  the  value  of  the  subscript  to  the  kth  element  in  the  molecular 
formula  of  species  i;  N^  is  the  number  of  moles  of  specie  i;  is  the 
fixed  number  of  moles  of  the  kth  element  in  the  system  (determined  from  the 
initial  conditions);  m  is  the  number  of  elements;  and  n  is  the  number  of 
species  (Smith  and  Missen,  1982,  p.  15). 

When  ions  are  present,  we  must  add  a  charge  balance  equation 
(conservation  of  electrons)  to  the  elemental  abundance  equations.  This 
equation  can  be  written  as 

n 

I  A^.N.  -  =  0;  k  =  m  +  1  (3.2) 

i=l 


where  the  coefficients  A^,  .,  i  =1, 

nrr  I  ,  1 

specie.  For  molecular  species,  A^^  . 


..,  n,  are  the  charge  of  the  ith 
=  0;  for  an  ion  with  a  +2  charge. 
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*nH-l  i  ~  ®  charge,  .  =  -1 ,  etc.  Since  the 
overall  charge  of  the  mixture  is  zero,  =  0.  The  complete  set  of 
constraint  equations  are 


n 

I  =  0;  k  =  1,  2,  ...,  m  +  1  (3.3) 


These  equations  can  also  be  written  in  matrix  form  as 

[A]  [M]  =  [B]  (3.4) 

where  [A]  is  the  formula  matrix,  [N]  is  the  species-abundance  vector,  and  [BJ 
is  the  element  abundance  vector  (Smith  and  Missen,  1982,  p.  16).  The  matrices 
may  be  written  in  detail  as 


A  A 

fll2  ... 

A21 

A22 

[A]  = 

• 

1 

^^m+1,1  ^nH-l,n^ 

j 

"2 

»2 

[N]  = 

(3.6)  and  [B] 

• 

• 

K. 

The  formula  matrix  [A]  must  be  checked  to  ensure  that  it  does  not  contain 
any  linearly  dependent  rows.  If  any  of  the  rows  are  linearly  dependent,  a 
singular  coefficient  matrix  will  result  (Smith  and  Missen,  1982,  p.  213).  The 
offending  rows  must  be  identified  and  eliminated.  One  technique  is  to  use 
Gauss-Jordan  reduction  to  reduce  the  left  side  of  [A]  to  unit  matrix  form 
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(Smith  and  Missen,  1982,  p.  24).  Linearly  dependent  rows  will  appear  as  a 
zero  vector  and  may  be  removed  from  the  system. 


3.2  Deriving  the  Lagrangian  Function 

We  are  now  ready  to  derive  the  equations  used  in  the  direct  minimization 
technique.  The  extremum  function  G  must  be  minimized,  subject  to  the  side 
conditions  of  equation  3.3.  This  type  of  constrained  optimization  problem  can 
be  solved  using  the  method  of  Lagrange  multipliers  (for  a  brief  review  of  the 
method,  see  van  Zeggeren  and  Storey,  1970,  p.  158).  The  Lagrangian  function  L 
is  defined  as: 

fiH-l  n 

L  =  G  +  I  ».  (  I  A.  .N.  -  8.  )  (3.8) 

k=1  i=l  ^ 

where  are  Lagrangian  multipliers.  Substituting  in  equation  2.19  for  G, 
we  get: 

n  m+1  n 

L  =  I  m.N.  +  I  ».  (  I  A.  .N.  -  B.  )  (3.9) 

i=1  ^  ^  k=l  i=l  ^  ^ 


Taking  the  derivative  and  setting  equal  to  zero: 
n  m+l 

dL  =  Z  OL/aN.)  dN.  +  I  (3L/air^)  d*.  =  0 
i=l  ^  ^  k=l 

Then,  since  ,  N^,  ...;  and  ...;  are  independent,  the 

corresponding  partial  derivatives  must  also  be  equal  to  zero 


OL/aN.),^^ 


m+1 

=  **i  ^  \i*k  = 


i  =  1 . n 


. . . ,  m  +  1 


(3.10) 


(3.11) 


(3.12) 
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Dividing  equation  3.11  by  RT  to  and  defining  =  w^/RT,  we  then 
have  the  set  of  equations 

m+1 

p./RT  +  =0;  1=1 . n  (3.13) 

and  the  original  conservation  equations 
n 

I  A.  .N.  -  8^  =  0;  k  =  1 . m  +  1  (3.14) 

i=l 

a  total  of  n  +  m  +  1  equations.  For  equilibrium  conditions  at  a  specified 
temperature  and  pressure  (T,P),  these  equations  can  be  solved  for  the  n  +  m  + 

1  unknowns  N.;  1  =1,  2 . n  and  k  =  1,  2 . m+1. 

For  equilibrium  conditions  at  a  specified  enthalpy  and  pressure  (aH,P), 
the  final  mixture  temperature  becomes  an  unknown  and  we  must  add  another 
equation  to  the  above  set.  This  equation  is  called  the  energy  equation  and 
can  be  written  as 

Hp(T,P,N^,  ...)p  =  H^d.P.N^.  ...)^  +  aH  (3.15) 

where  Is  the  total  enthalpy  of  the  reactants,  Hp  Is  the  total  enthalpy 
of  the  equilibrium  products  and  AH  is  some  specified  enthalpy  change.  The 
temperature  is  solved  by  Iteration  like  any  of  the  other  unknowns.  In  many 
cases  the  enthalpy  change  may  be  due  to  heat  transfer;  for  the  adiabatic  case, 
AH  =  0. 

The  set  of  equations  3.13,  3.14,  and  3.15  contain  several  parameters  that 
we  must  derive  expressions  for.  These  Include  the  chemical  potential  of  each 
species  (p^.),  the  product  and  reactant  enthalpy  (Hp,  H^),  and  the 
number  of  moles  of  each  element  In  the  system  (Bj^). 
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3.3  Enthalpy 

The  total  enthalpy  of  a  mixture  (whether  it  be  a  mixture  of  products  or  a 
mixture  of  reactants)  may  be  calculated  from 

n 

H  =  I  H.N.  (3.16) 

i=l  ^  ^ 

where  is  the  partial  molal  (or  molar)  enthalpy  of  species  i.  The  partial 
molar  enthalpies  of  the  species  can  be  calculated  from  one  of  the  following 
equations  depending  on  the  convention  used  to  specify  the  activity 
coefficients  of  the  species  (Denbigh,  1981,  p.  128  and  p.  279).  For  the 
convention  -y .  1  as  -*  1 ,  then 

H.  =  h^  -  RT^(alnY^/3T)  (3.17) 

where  h.  is  the  enthalpy  per  mole  of  the  pure  species  i.  For  the  convention 
-»  1  as  ■»  0  and  the  convention  y^  1  as  m^  -»  0,  then 

=  h’  -  RT^{3lnY/3T)  (3.18) 

O 

where  H.  is  the  partial  molar  enthalpy  of  species  i  at  infinite 
dilution.  Generally,  the  partial  derivatives  in  the  above  equations  are 
calculated  numerically  by  perturbing  the  appropriate  function  with  respect  to 
temperature.  Of  course,  for  an  ideal  solution,  the  activity  coefficients  are 
always  equal  to  one,  and  the  second  term  in  each  of  the  above  equations  will 
disappear. 

o 

The  temperature  dependence  of  h^  and  H.  can  be  calculated  from 
specific  heat  data.  Experimental  data  for  the  specific  heat  of  a  substance 
can  be  curve-fit  to  a  polynomial  of  the  form 
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(3.19) 


Si  ■  ^ 

where  Cp^  is  heat  capacity  per  unit  mole  of  the  pure  component  i,  or 

Si  “  ^i  ‘’i^  *  ‘^i^^  (3.20) 

O 

where  Cp^  is  the  partial  molar  heat  capacity  of  species  i  at  infinite 

o 

dilution.  Using  either  the  equation  dh^  =  Cp^dT,  or  the  equation  dH. 

o 

=  Cp.dT,  the  polynomial  can  be  integrated  to  get  the  enthalpy  of  the 
species  as  a  function  of  temperature: 

h,  =  J  C  .dT  +  constant 

1  pi 

o  o 

H.  =  j'  C  .dT  +  constant 
1  pi 

or 

h.  =  a.T  +  b^T^/2  +  c^T^/3  +  d^  (3.21) 

h’  =  a.T  +  b^T^/2  +  C.T^/3  +  d^  (3.22) 

The  integration  constant  d^  can  be  calculated  if  the  enthalpy  is  known  at 
any  temperature  within  the  range  of  applicability  of  the  equation. 

3.4  Chemical  Potential 

From  Denbigh  (1981)  the  chemical  potential  for  a  single  perfect  gas  is 
defined  as  follows: 

p(T,P)  =  n"(T)  +  RTln  P  (3.23) 

where  w"(T)  is  the  standard  chemical  potential  (or  Gibbs  free  energy  per 
mole)  at  temperature  T  and  one  atmosphere  pressure;  P  is  the  pressure  in 
atmospheres,  and  R  is  the  universal  gas  constant. 
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For  a  single  Imperfect  gas 


ti(T,P)  =  n  (T)  +  RTln  f, 

f  =  yP,  (3.24) 

Y  -*  1  as  P  -»  0 

where  f  is  the  fugacity  of  the  gas,  y  the  fugacity  coefficient,  and 
p®(T)  is  the  chemical  potential  of  the  gas  at  unit  fugacity.  Since  y  is 
very  nearly  1  at  P  =  1 ,  the  numerical  value  for  ii“(T)  is  essentially  the 
same  as  in  equation  3.23.  Note  that  for  a  perfect  gas,  y  =  1  and  equation 
3.24  reduces  to  equation  3.23. 

For  a  perfect  gas  mixture 

P^(T.P.X^)  =  vJ(T)  +  RTln  P.  (3.25) 

where  u^(T)  is  the  standard  chemical  potential  of  gas  i  in  its  pure 
state  at  atmospheric  pressure  (exactly  the  same  as  u^(T)  in  equation 
3.23).  is  the  mole  fraction  of  component  i  in  the  gaseous  mixture 
defined  by 

ngas 

X.  =  N./  I  N.  (3.26) 

1  i=l  ^ 


and  P.  is  the  partial  pressure  of  gaseous  species  i  defined  by  P.  =  X.P. 
For  an  imperfect  gas  mixture 


u^d.P.X.)  =  w.(T)  +  RTln  f^ 


f,  = 


(3.27) 


where  f.  is  the  fugacity  of  species  i,  y.  is  the  fugacity  coefficient  of 
species  i  in  the  gaseous  mixture,  and  w^(T)  has  the  same  value  as  in 


17 


equation  3.24.  Once  again,  for  an  Ideal  gaseous  mixture  =  i  and 
equation  3.27  reduces  to  equation  3.26. 

The  chemical  potential  for  condensed  phase  solutions  is  defined 
similarly.  When  no  distinction  is  made  between  the  solvent  and  solutes  in  a 
solution,  the  Raoult  convention  is  commonly  used  for  all  species  in  that 
solution.  When  a  distinction  is  made  between  solvent  and  solutes  (e.g., 
aqueous  solutions),  the  Henry  convention  is  commonly  used  for  the  solute 
species  and  the  Raoult  convention  for  the  solvent  species  (Smith  and  Missen, 
1982,  p.  52). 

For  a  Raoults'  law  ideal  solution 

P^(T,P,X^)  =  u*(T,P)  +  RTIn  (3.28) 

■k 

Where  n^(T,P)  is  the  chemical  potential  of  pure  species  i  at  the 
temperature  and  pressure  of  the  solution. 

For  a  solute  following  the  Henrys*  law  convention 

y^(T,P.X^)  =  V*(T,P)  +  RTIn  X^  (3.29) 

A 

where  v^(T,P)  is  the  chemical  potential  of  solute  species  i  in  an 
infinitely  dilute  solution  and  X^  is  te  mole  fraction  of  the  ith  liquid 
species  in  the  particular  liquid  phase. 

For  a  nonideal  solution 

W^(T.P,X.)  =  y*(T,P)  +  RTIn  a^  (3.30) 

where  a^  is  the  activity  of  species  i  in  the  solution.  For  non-electrolytic 
solutions  the  activity  can  be  written  as 

a^  =  Y^X.  (3.31) 
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where  the  activity  coefficient  of  species  i  in  the  solution  which 

accounts  for  any  deviation  from  ideality.  Note  that  for  an  ideal  solution, 

=  1 .  For  a  species  which  approaches  Raoults'  law,  y^  1  as 
-*  1.  For  a  solute  species  which  exhibits  Henrys'  law  behavior,  y^  ■*  1 
as  X.  -»  0. 

In  some  cases  the  amount  of  solute  is  specified  on  the  molality  scale. 

Then 

w^(T,P,m^)  =  p*(T,P)  +  RTln  a^ 

(3.32) 

Y^  -♦  1  as  m.  -♦  0 

where  a^  =  y^*"^  and  m^  =  1000  N^/(M^N^).  is  the  molecular  weight  of  the 

■k 

solvent  and  is  the  number  of  moles  of  the  solvent.  >i^(T,P)  is  the  chemical 
potential  of  the  solute  in  a  hypothetical  ideal  solution  of  unit  molality  at 
the  same  temperature  and  pressure  as  the  solution  under  discussion  (Denbigh, 
1981,  p.  276). 

The  chemical  potential  of  a  single  species  phase  such  as  a  pure  solid  is 
written  as 

»«^(T,P)  =  u*(T,P)  (3.33) 

* 

where  n^(T,P)  is  the  chemical  potential  of  the  pure  species  i  at  the 
temperature  and  pressure  of  the  phase.  All  of  the  above  specifications  for 
chemical  potential  contain  a  standard  reference  potential  for  that  species. 

It  is  necessary  to  assign  a  numerical  value  to  that  term.  According  to  van 
Zeggeren  and  Storey,  (1970),  p.  31,  we  may  set 

gJ(T)  =  aGj.(T)  (3.34) 

-  19  - 


where  &G^^(T)  is  the  standard  Gibbs  free  energy  of  formation  of  the 
pure  species  (or  ideal  one  molal  solution,  whichever  the  case  may  be)  at 
temperature  T.  Similarly,  for  reference  potentials  that  are  a  function  of 
pressure 

g*(T,P)  =  AG*.(T.P)  (3.35) 

* 

where  &Gp.(T)  is  the  standard  Gibbs  free  energy  of  formation  of  species 
i  at  temperature  T  and  pressure  P.  These  quantities  can  be  found  tabulated 
for  a  wide  range  of  species  in  the  literature  e.g.,  Wagman  et  al.  (1968)  and 
JANAF  (1971).  If  it  is  only  available  at  one  temperature,  values  at  other 
temperatures  can  be  calculated  from  enthalpy  data  using  a  form  of  the  van't 
Hoff  equation  (Denbigh.  1981,  p.  143): 

[a(-aGf ./T)/aT]p  =  (3.36) 

where  the  enthalpy  of  formation  represents  the  enthalpy  change 
occurring  when  species  i  is  formed  from  the  elements  comprising  it.  For  a 
solvent,  a  gas  or  a  Raoults'  law  solute,  AHj^^  may  be  calculated  from 
the  enthalpy  per  mole  of  the  pure  component,  h.  (equation  3.21),  and 
enthalpy  per  mole  of  the  elements  comprising  it 

Ah“  =  h.  -  Z  v.h.  (3.37a) 

*  '  *  p  V  U 

where  hj  and  vj  are  the  enthalpy  and  the  stoichiometric  coefficient  of 
the  jth  element  in  the  formation  reaction.  For  a  Henry's  law  solute, 

O 

AHp.  may  be  calculated  from  the  partial  molar  enthalpy  of  the  species 

o 

at  infinite  dilution,  H^  (equation  3.22),  and  the  enthalpy  per  mole  of 
the  reference  elements: 
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(3.37b) 


Assuming  that  the  enthalpies  of  the  elements,  h^,  are  fitted  to  the 
same  type  of  function  as  the  species  formed  (given  by  equations  3.21  and 
3.22),  we  have 

=  A^T  +  B^T^/2  +  C.T^/3  +  0^  (3.38) 

where  A.  =  a^  -  X  ®i  -  ^1  “  ^  I'’t®9'’3ting  equation  3.36, 

-agJ^/T  =  J  (A^/T  +  B^/2  +  C^T/3  +  D./T^)  dT  +  constant 
Finally, 

=  A^T  ln(T)  -  B.T^/2  -  C.T^/6  +  0^  +  E^T  (3.39) 

Since  most  compilations  list  the  Gibbs  free  energy  of  formation  at  a 
pressure  of  one  atmosphere  only,  we  similarly  can  calculate  values  at  other 
pressures  using  (Denbigh,  1981,  p.  300): 

[3(Ag“./T)/3P].j.  =  AV*.  (3.40) 

9 

where  AV^^  is  the  standard  vaolume  change  for  the  formation  of  species 

o  o 

i,  analogous  to  AH^^  and 

3.5  Specifying  Amounts  of  Reactants 

The  relative  amounts  of  reactants  present  can  be  specified  by  a  parameter 
called  the  mixture  fraction  (Spalding,  1979,  ch.  6),  usually  denoted  by  the 
symbol  f.  Suppose  we  divide  the  components  of  the  reactants  into  two 
portions:  the  "O"  (zero)  fluid  and  the  "®"  (infinity)  fluid.  The  fluid 
resulting  from  the  mixing  of  the  0  fluid  and  the  <»  fluid  will  be  called  the 
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"p"  (product)  fluid.  For  any  conserved  property  <t>,  the  mixture  fraction  is 
defined  as 


f  =  (♦P  -  -  ♦“)  (3.41) 

Usually,  the  fluid  state  at  »  and  0  is  known,  and  we  wish  to  determine 
at  a  given  mixture  fraction.  Rewriting  equation  3.41,  we  get 

4>P  =  f4>°  +  (1  -  f)  (3.42) 

Therefore,  for  f  =  1 ,  «|)P  =  ,  and  for  f  =  0,  =  «j>“.  If 

we  exclude  nuclear  reactions,  then  Yj^,  the  mass  fraction  of  element  k  in  a 

mixture  of  compounds,  is  a  conserved  property  (Kays  and  Crawford,  1980,  p. 

344).  Then  we  can  write 

=  fY°  +  (1  -  f)  k  =  1,  2 . m  (3.43) 

where 


Y 

Y 


0 

k 


k 


k  =  1,  2,  ...,  m  (3.44) 

k  =  1,  2 . m  (3.45) 


0  (D 

Y. ,  Y^ ,  are  the  mass  fraction  of  species  i  in  the  0  fluid  and  « 
fluid,  respectively;  and  H.  are  the  molecular  weights  of  element  k  and 
species  i  respectively;  and  Aj^.  are  the  number  of  moles  of  element  k  in 
species  i  (defined  previously). 

The  number  of  moles  of  element  k  in  the  system,  Bj^,  can  be  written  as 


“k ' 


Using  equation  3.43,  we  get 


k  =  1 ,  2,  . . . ,  m 


(3.46) 
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\  '  +  (1  -  n 


k  =  1 ,  2,  . . . ,  m 


“k '  f  \iV"i  *<'-'> 


k  =  1.  2 . tn  (3.47) 


Equation  3.47  can  be  used  to  calculate  the  components  of  the  B  vector  when  the 
mixture  fraction  f,  the  mass  fractions  of  all  species  in  the  0  fluid  state 

o  o 

Y.  and  in  the  ®  fluid  state  Y.  are  specified.  For  the  B  vector, 
the  computer  code  developed  in  later  section  take  the  mass  of  the  mixture  as 
one  kg  such  that  Bj^  can  then  be  interpreted  as  the  number  of  moles  of  kth 
element  in  a  unit  kg  mass  of  mixture. 
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4.  NUMERICAL  METHOD 


4.1  Newton-Raphson  Technique 

The  system  of  nonlinear  equations  set  up  by  the  direct  minimization 
technique  in  Section  3  can  be  solved  using  Newton-Raphson  iteration  (Gerald 
and  Wheatley,  1984,  pp.  135-137).  The  set  of  r  =  n  +  m  +  1  equations  to  be 
solved  can  be  written  in  the  following  form: 

(X-l .  ^2,  . . . ,  X^)  =  0 
^2  X2 . X^)  =  0 

(4.1) 


fp  (X^.  X2 . X^)  =  0 


where  X  ,  X2,  ...,  X^  are  the  roots  of  the  equations.  The  relationship 
between  any  f^  (X^  +  4X^ ,  X  ^  +  6X^,  ...)  and  f.  (X^,  X2,  ...)  where 
iX^ ,  4X2,  ...,  4X^  denote  arbitrary  increments  in  X.|,  X2,  ...,  X^  is  given  by 
Taylor's  Expansion  (Callen,  1985,  p.  474): 

f.  (X^  +  AX^,  X2  +  4X2,  ...) 

=  fi  (X,.  X2,  ...)  > 

+  (af^/ax^)  4X^  +  (3f./3X2)  4X2  +  ...]  (4.2) 

+  (1/2)  [(a^f./3X^)  4X^  +  ...) 

+  •••;  i  =  1.  2 . r 

If  we  define  4X.  =  X^  '  ^i’  4nd  X.  are  the  true  and  guessed  roots, 

respectively,  so  that 


f,  (X, 


+  4X^  ,  X2  +  4X2, 


)  =  (X^.  X2,  ...)  =  0; 


i  =  1.  2, 


(4.3) 
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we  could  solve  the  above  set  of  equations  for  the  4X.  and  calculate  the 
roots  from  X^  =  X^  +  4X..  In  practice,  if  the  initial  estimates 
X^  are  near  enough  to  the  solution,  we  can  truncate  the  series  after  the 
first  order  terms  and  get  an  approximation  X  to  the  final  solution  X.: 

f^  (x'[^\  x^\  ...)  =  0  =  f^  (x™,  X™,  ...)  +  [af/ax^)  ax^  +  ...] 

or 

-  f^  (x']’,  x^.  ...)  =  (af./ax^)  ax^  (af./ax^)  ax^  + 

1=1,2 . r  (4.4) 

where  x"^^  =  X™  +  ax.  This  set  of  linear  equations  is  solved  using 
Gaussian  elimination  (Gerald  and  Wheatley,  1984,  p.  91)  for  the  unknowns 
ax^.  The  new  approximations  are  calculated  from  X*"^^  =  X™  +  aX  and 
the  procedure  is  repeated  until  or  f  (X™,  X^, 

where  and  are  some  small  numbers  {l.E-6  and  5.E-5  respectively, 
in  this  case) . 

In  order  to  solve  the  above  set  of  equations  we  must  calculate  the 
partial  derivatives  appearing  on  the  right  sides.  Rather  than  derive  an 
expression  for  the  partial  derivatives  by  hand,  it  is  easier  to  approximate 
the  partials  numerically.  This  is  done  by  recalculating  the  function  with  a 
small  perturbation  to  each  of  the  variables  in  turn  (Gerald  and  Wheatley, 

1984,  p.  137); 


af/ax^  =  [f^ 

(X^  +  4.  x^. 

...)  -  f^  (X^,  x^,  .. 

.)]/4 

af^/aXj  =  [f^ 

(x^.  x^  +  6, 

...)  -  f^  (X^,  x^,  .. 

.)]/< 

(4.5) 

af^/ax^  =  [f^^  (X^ . x^  +  a)  -  f^  (x^ . x^)]/a 
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Thus,  the  solution  to  the  set  of  nonlinear  equations  Is  approximated  by 
solving  the  linear  set  of  equations.  Each  successive  approximation  generates 
a  set  of  solutions  closer  and  closer  to  the  final  solution. 

4.2  Relaxation  Factor 

Certain  modifications  to  this  technique  are  necessary  in  order  to  improve 
the  stability  of  the  iteration.  A  step-size  parameter  (or  relaxation  factor) 

T  can  be  introduced  as  follows: 

xiiH-l  ^  (4.6) 


where  0  <  t  <  1.  Poor  initial  guesses  can  lead  to  a  large  negative 
correction  which  could  cause  the  number  of  moles  of  a  species  to  become 
negative.  We  choose  the  relaxation  factor  to  ensure  that  the  mole  numbers  of 
all  the  species  remain  positive  throughout  the  iteration  process 
(non-negativity  constraint).  We  start  the  process  with  t  =  1  and  find  the 
maximum  value  for  t  less  than  one  that  will  ensure  that  all  of  the  new  mole 
number  estimates  will  be  positive.  The  algorithm  used  is 


Start  with  T  =  1 

For  i  =  1 ,  2 . n 

If  <  0,  Then 

(4.7) 

Tmax  =  -  c)/ax^) 


T 


'max 


Next  i 


where  c  is  some  small  number  (0.01  in  this  case).  Note  that  there  is  no 
such  non-negativity  constraint  on  the  Lagrangian  multipliers.  The  Euclidean 
norm  of  the  error  vector  is  defined  as  (Gerald  and  Wheatley,  1984,  p.  115): 
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(4.8) 


II  f"*  IL  =  [  X  U  (X?,  Xp* 

e  i=l  '  ' 

As  long  as  the  Newton-Raphson  corrections  decrease  the  norm  of  the  error 
vector  (i.e.,  the  iteration  is  converging)  t  is  set  equal  to  If 

for  any  iteration  m  +  1,  the  norm  is  greater  than  the  norm  calculated  at  the 
previous  step,  i.e. 

II  f"^^  II  >  II  f"*  II  (4.9) 

the  relaxation  factor  is  reduced  by  setting  ^  equal  to  and  the 

norm  at  the  m  +  1  step  is  recalculated  (see  Prausnitz  et  al.,  1980,  p.  116). 
This  reduction  of  t  is  repeated  until  the  norm  decreases  (in  which  case  the 
iteration  proceeds),  or  until  i  becomes  too  small,  say  0.01,  (in  which  case 
T  is  just  set  equal  to  and  the  iteration  is  allowed  to  take  a  step 
which  increases  the  norm  of  the  error  vector). 


4.3  Species  Present  in  Small  Amounts 

When  a  species  in  a  multispecies  phase  is  present  in  a  very  small  amount, 
it  forces  the  step-size  parameter  to  become  very  small,  slowing  the 
convergence  of  the  algorithm.  The  species  must  then  be  removed  from  the  main 
calculation.  After  convergence  is  achieved  with  this  species  removed,  the 
amount  of  the  species  present  may  be  calculated  form  (Smith  and  Missen,  1982, 
pp.  217-218): 


X 


i 


*  nn-l 
exp  (-  n^/RT  - 


(4.10) 


This  equation,  of  course,  is  only  valid  if  the  species  is  indeed  present  in  a 
very  small  amount  at  equilibrium.  Note  that  during  the  course  of  iteration,  a 
species  may  be  erroneously  removed  (especially  if  poor  initial  estimates  are 
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given).  If  so,  equation  4.10  provides  a  means  of  checking  If  any  other 
species  should  have  been  Included  In  the  main  calculation. 

An  equation  similar  to  that  given  above  may  be  used  to  determine  If  a 
single  species  phase  (such  as  a  pure  solid  species)  should  be  Included  In  the 
calculation.  Equilibrium  may  at  first  be  calculated  with  only  raultlspecles 
phases  present,  then  the  following  test  may  be  made  to  determine  If  the 
conditions  are  favorable  for  the  formation  of  a  single  species  phase  (Smith 
and  Hissen,  1982,  p.  59): 

* 

p./RT  +  <  0  (4.11) 

If  the  quantity  given  above  is  less  than  zero,  the  single  species  phase  in 
question  should  be  added  to  the  system  and  equilibrium  recalculated.  If  the 
quantity  is  exactly  zero,  the  phase  Is  at  incipient  formation.  If  the 
quantity  Is  greater  than  zero,  the  phase  should  not  be  Included.  A  similar 
test  may  be  used  for  multispecies  phases,  as  well  (Smith  and  Missen,  1982,  p. 
59). 
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5.  COMPUTER  PROGRAM 


The  equations  developed  In  the  previous  sections  have  been  implemented  in 
a  computer  program  called  CEC-NMS  (Complex  Equilibrium  Calculations  of 
Nonideal,  Multiphase  Systems).  For  a  listing  of  the  program,  see  the 
appendix.  The  program  is  written  in  Fortran  77  and  is  approximately  1300 
lines  long,  including  a  generous  number  of  comment  cards.  CEC-NMS  consists  of 
a  main  program  unit  and  9  primary  subroutines.  Additional  subroutines  to 
calculate  the  activity  coefficients  for  nonideal  solutions  may  also  be 
present,  depending  on  the  particular  system. 

5.1  Main 

The  main  program  unit  performs  several  functions: 

1.  Declare  and  initialize  variables. 

2.  Set  up  chemical  system  (species  symbols,  molecular  weight,  species 
coefficient  matrix,  etc.). 

3.  Open  output  file  and  read  in  data  for  calculation  point. 

4.  Call  subroutine  GUESS  to  provide  initial  estimates. 

5.  Calculate  reactant  enthalpy. 

6.  Call  subroutine  NLSYST,  the  Newton-Raphson  equation  solver. 

7.  Print  out  error  messages,  if  any. 

8.  Call  subroutines  PREOUT  and  OUTPUT  to  print  out  results. 

9.  Return  to  step  3  and  read  in  data  for  the  next  calculation  point. 
Stop  program  execution  after  last  data  point. 

5.2  Gcal 

Subroutine  GCAL  calculates  the  standard  Gibbs  potential  of  each  species 
in  the  system  at  a  given  temperature.  In  the  present  form,  equation  3.39  is 
used.  The  user  must  supply  the  coefficients  A^,  B^,  etc. 
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5.3  Heal 

Subroutine  HCAL  calculates  the  partial  nwlar  enthalpy  of  the  species  in 
solution  (H^  in  equation  3.16)  using  equations  3.17-18  and  3.21-22.  The 
user  must  supply  the  coefficients  a^,  b^,  etc.  Later  versions  of  the 
program  combine  subroutine  GCAL  and  subroutine  HCAL  into  one  subroutine  called 
HSGCAL. 

5.4  Preout 

Subroutine  PREOUT  calculates  the  number  of  moles  of  any  species  which 
have  been  removed  from  the  main  calculation  using  equation  4.10.  Also 
calculated  are  various  mass  and  mole  fractions,  the  mixture  density,  and  the 
mixture  void  fraction  (volume  fraction  of  gas  phase  present).  The  mixture 

3 

density,  p  [kg/m  ],  is  defined  as 

P  =  1/V  (5.1) 

3 

where  V  is  the  specific  volume  of  the  mixture  [m  /kg].  The  specific  volume 
of  a  mixture  of  L  phases  may  be  expressed  as 

L 

V  =  I  Y.V.  (5.2) 

j=l  J  ^ 

or 

L 

V  =  X  Y,/p.  (5.3) 

j=l  J  ^ 

where  V^,  pj,  and  Yj  are  the  specific  volume,  density,  and  mass 

fraction  of  phase  j  in  the  total  mixture,  respectively.  Substituting  equation 

5.3  into  5.1 ,  we  have 

L 

P  =  l/(  X  Y/p.)  (5.4) 

j=1  ^  J 
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The  void  fraction  of  the  mixture  is  defined  as 


a 


(5.5) 


or 


ngas  n 

a  =  (  X  Y/p.)/(  I  Y./p.)  (5.6) 

j=i  ^  ^  j=l  J 


5.5  Output 

Subroutine  OUTPUT  sends  detailed  results  for  a  calculation  point  to  the 
line  printer  and  abridged  results  to  a  file.  The  file  is  later  used  to  plot 
the  results. 

5.6  Guess 

Subroutine  GUESS  provides  initial  estimates  needed  to  begin  the  iteration 
at  each  calculation  point.  The  user  must  provide  these  estimates  for  every 
species.  Fortunately,  CEC-NMS  has  converged  quite  well  for  all  test  systems 
with  rather  arbitrary  initial  guesses. 

5.7  Nlsvst 

Subroutine  NLSYST  solves  the  system  of  nonlinear  equations  by  a 
step-limited  Newton-Raphson  technique.  The  subroutine  performs  the  following 
functions: 

1.  Declare  and  initialize  variables. 

2.  Begin  calculations  with  only  multispecie  phase  present.  Single 
species  phases  will  be  added  one  at  a  time  after  convergence  is 
achieved  with  only  multispecies  phases  present. 

3.  Call  subroutine  FCN  to  calculate  function  values. 

4.  Remove  any  specie  or  phase  that  is  present  in  a  very  small  amount 
( i .e. ,  less  than  lO”^®) . 
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5.  Remove  any  specie  that  Is  causing  convergence  problems  (i.e.,  if 
convergence  is  not  achieved  within  20  iterations  since  the  last 


species  was  removed,  remove  the  species  in  the  system  that  is 
present  in  the  1(  ^t  amount). 

6.  Cueck  if  function  values  meet  the  F-tolerance. 

7.  Calculate  the  partial  derivative  matrix. 

8.  Call  subroutine  SCALE  to  scale  the  matrix. 

9.  Call  subroutine  ELIM  to  solve  the  matrix. 

10.  Calculate  the  maximum  relaxation  factor  that  will  ensure  that  the 
mole  numbers  of  all  species  remain  positive. 

11.  Reduce  the  relaxation  factor  if  such  a  reduction  will  reduce  the 
Euclidean  norm  of  the  vector. 

12.  Calculate  new  estimates  for  the  unknowns. 

13.  Check  if  X-tolerance  is  met. 

14.  After  convergence  is  achieved  with  only  multispecies  phases  present, 
check  if  any  of  the  single  species  phases  will  reduce  the  Gibbs  free 
energy  of  the  system  using  equation  4.11.  Recalculate  the 
equilibrium  with  each  new  phase  added  in  turn. 


5.8  Elim 

Subroutine  ELIM  solves  the  partial  derivative  matrix  generated  by  NLSYST 
(the  set  of  equations  4.4)  using  Gaussian  elimination. 

5.9  Scales 

This  subroutine  scales  the  values  of  the  partial  derivative  matrix  so 
that  the  largest  element  in  each  row  is  unity. 
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5.10  Fen 


Subroutine  FCN  calculates  the  function  values  (the  in  equations  4.4) 
of  the  system  of  nonlinear  equations.  The  following  calculations  are 
performed: 

1.  Shift  the  X-vector  if  any  species  have  been  removed. 

2.  Calculate  the  enthalpy  of  the  mixture  by  calling  subroutine  HCAL. 

3.  Calculate  the  standard  Gibbs  free  energy  of  the  species  by  calling 
subroutine  GCAL. 

4.  Calculate  the  mole  fractions  of  each  species. 

5.  Calculate  the  chemical  potential  of  each  species  using  equations 
3.25,  3.27,  3.28,  3.30,  or  3.32. 

6.  Calculate  the  functions  using  equations  3.13-15. 

7.  Shift  F-vector  for  any  removed  species. 

5.11  Activ 

Optional  user  supplied  subroutines  ACTIVl,  ACTIV2,  etc.  are  used  to 
calculate  the  activity  coefficients  of  multispecies  phase  1,  phase  2,  etc. 
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6.  APPLICATION  TO  AN  ELECTROLYTIC  SYSTEM 


In  the  course  of  study  of  the  structure  of  chemically  reacting  jets  it 
became  necessary  to  compute  the  equilibrium  compositions  for  a  mixture  of 
aqueous  ammonia  (NH^)  solution  and  hydrogen  chloride  gas  (HCl)  (Chan  et  al... 
1987).  These  reactants  form  a  multicomponent  electrolytic  solution  containing 
ions  and  molecules  dissolved  in  water. 


The  initial  compositions  of  the  reactants  are  specified  by  the  mixture 

fraction  f  (see  equation  3.41)  and  the  mass  fractions  of  species  in  0  and  « 

0  ^ 

fluid  streams,  Y.  and  Y..  The  0  and  *  fluids  were  defined  to  be 

0  OB 

pure  HCl(g),  (Y.  =  1),  and  aqueous  ammonia  solution 

respectively.  The  initial  concentration  of  the  aqueous  ammonia  solution 
(Y“H3(aq))  varies  form  0-30%  by  mass.  The  relative  amount  of  HCl  to 
ammonia  solution  ranged  from  0-100%  (i.e.,  the  mixture  fraction  ranged  from  0 
to  1). 


The  system  under  consideration  was  assumed  to  contain  the  following 
species  in  the  liquid  phase:  NH^,  Cl”,  OH*  and 

^2^(L)'  chemical  potential  of  the  solutes  were  expressed  using 
equation  3.32.  the  chemical  potential  of  water  (the  solvent)  was  expressed 
using  equation  3.30.  The  gas  phase  (when  present)  was  assumed  to  consist  of 
^2*^(g)’  ^^3(g)'  chemical  potentials  of  the  gaseous 

species  were  expressed  using  equation  3.27.  A  solid  phase  consisting  of 
NH^Cl^c)  appears  at  high  concentrations  of  NH^  and  Cl”.  The 
chemical  potential  of  the  pure  solid  phase  was  expressed  using  equation  3.33. 


6.1  Species  Coefficient  Matrix 

The  formula  matrix  [A]  of  the  system  is  written  as 
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where 
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(E) 

(F) 
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(H) 

(I) 

(J) 
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0 
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0 
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1 

-1 

-1 

0 

0 

0 

0 

0 

0 

A  = 

0  = 

0H~ 

G 

=  H20 

(g) 

j 

=  nh^ 

"(c) 

B  =  NH^ 

E  = 

"“scaq) 

H 

=  "«3(9) 

C  =  Cl' 

F  = 

H2O 

(L) 

I 

=  HCl 

(g) 

and 

Cl  =  chlorine 
H  =  hydrogen 
N  =  nitrogen 
0  =  oxygen 
e  =  electron 

Performing  Gauss-Jorden  reduction  we  obtain 


(A)  (B)  (C)  (0)  (E)  (F)  (G)  (H)  (I)  (J) 


H 

N 

Cl 

0 

e 


1 
0 
0 
0 
10 


0 

1 

0 

0 

0 


0 

0 

1 

0 
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0 

0 

0 

1 
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-1 

1 

0 
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0 

0 
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1 

0 

0 

1 

1 


-1 

1 

0 

0 
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1 

0 

1 

0 
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0 

1 

1 

0 
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We  see  that  row  4  and  row  5  are,  in  fact,  the  same  equation.  Therefore,  we 
eliminate  row  5,  the  charge  balance  equation,  as  it  is  redundant.  The  final 
revised  formula  matrix  is 
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(A) 

(B) 

(C) 
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(H) 

(I) 

(J) 

Cl 
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-1 
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1 

0 

0 

1 

0 

0 

1 

0 

1 

0 

vO 

0 

0 

1 

0 

1 

1 

0 

0 

0. 

Therefore,  for  constant  temperature  and  pressure  equilibrium,  there  are  a 
total  of  14  unknowns  -  the  mole  numbers  of  10  species  and  the  Lagrangian 
multipliers  for  4  constraint  equations.  Note  that  if  during  the  course  of 
computation,  any  species  is  eliminated  from  consideration,  we  must  check  to 
see  that  all  of  the  rows  in  the  matrix  are  still  independent. 


6.2  Calculation  of  Thermodynamic  Data 

The  thermodynamic  data  necessary  for  this  system  was  obtained  from  a 

o  o  o 

variety  of  sources.  C^,  AH^,  and  A6^  for  the  ions  at  298.15  K  is  from 

O 

Wagman  et  al.  (1968).  AHj.(T)  for  was  curve-fit  from  experimental 

data  found  in  Washburn  (1930).  The  polynomial  constants  for  the  specific  heat 
of  the  gases  were  obtained  from  Reynolds  and  Perkins  (1984).  Thermodynamic 
data  for  NH^Cl^^j  is  from  the  JANAF  Thermochemical  Tables  (1971).  The 
computer  program  used  to  calculate  the  constants  A.,  8.,  C.,  0.,  and  of 
equation  3.38  and  3.39  is  listed  in  Appendix  A. 


6.3  Activity  Coefficients  in  the  Liquid  Phase 

The  activity  coefficients  of  the  species  in  the  liquid  phase  were 
approximated  using  a  theoretical  model  developed  by  Edwards  et  al.  (1978). 
The  Edwards  model  is  based  on  formalism  set  forth  by  Pitzer  in  a  series  of 
papers  (1973,  1974).  Edwards,  however,  neglects  all  ternary  interaction 
parameters.  Their  justification  for  this  simplification  is  that  sufficient 
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ternary  experimental  data  is  rarely  available  to  fix  ternary  parameters. 
Therefore,  the  model  relies  instead  on  only  binary  parameters  to  describe  the 
interactions  between  the  species  in  solution.  Edwards  also  extended  the 
equations  to  weak  electrolytes  and  molecular  species  by  defining  parameters 
for  ion-molecule  and  molecule-molecule  interactions  as  well  as  the  Pitzer 
ion-ion  interactions. 

It  should  be  noted  that  Zemaitis  et  al.  (1986)  provides  an  excellent 
summary  of  the  most  recent  modeling  techniques  for  determining  the  activity 
coefficients  in  electrolytic  solutions.  Included  are  chapters  on  single  and 
multicomponent  strong  electrolytes,  strongly  complexing  compounds,  weak 
electrolytes  and  molecular  species.  The  equations  Edwards  used  for  the 
activity  coefficient  of  a  solute  are  (Zemaitis  et  al.,  1986,  pp.  503): 


In  y,  =  Z 


f^  +  2 


ns 


*  Z 


ns 

I 


ns 

1 


j=l  k=l 


'"j"’k®jk 


(6.1) 


where 


f^  =  -A^  [v'I/(l  +  l-Z/I)  +  (2/1.2)  In  (1  +  1.2/1)] 

♦ 

®ij  "  ®ij  ^  [1  -  (1  +  2^1)  exp  (-  2/1)] 

ejj  =  (B]j/4I^)  (-  1  +  (1  +  2/1  +  21)  exp  (-  2/1)] 

A.  =  natural  log  based  Oebye-Huckel  constant  (see  equation  6.3) 
9 

I  =  ionic  strength  of  the  solution 
ns  y 

=  0.5  I  Z‘m 
i=1  ’  ’ 

and  where 


-  37 


=  ionic  charge  of  species  i 


m.  =  molality  of  species  i 

=  binary  interaction  parameter 
®ij  “  binary  interaction  parameter 

ns  =  number  of  solutes  in  the  solution.  Note  that  this  does  not  include 
”2°(1)- 

The  activity  of  water  is  expressed  as  (Zemaitis  et  al.,  1986,  pp.  504); 


In  a^  =  {2A^I^^^/(1  +  1.2^I) 

(6.2) 

ns  ns  o  1  ns 

-II  [014  +  B-i  exp  (-  2/1)]  -  I  m.} 

i=l  j=i  1  J  u  i=l  ^ 


The  Oebye-Huckel  parameter  for  aqueous  electrolytic  systems  is  calculated 
from  (Chen  et  al.,  1982): 


=  -  61.44534  exp  [(T  -  T  )/T  ] 

9 

+  2.864468  [exp  (T  -  T*)/T*]^ 

4-  183.5379  In  (T/t“)  -  0.6820223  (T  -  t’) 

0.0007875695  (T^  -  (T*)^)  +  58.95788  (t’/T)  (6.3) 

O 

where  T  =  273.15  K.  The  estimated  values  for  the  binary  interaction 
parameters  used  in  this  study  along  with  the  source  are  shown  in  Table  1. 

Note  that,  as  suggested  by  Edwards,  like  charge  ionic  interactions  are 
ignored,  as  well  as  the  b1.  term  for  interactions  involving  molecules. 

*  J 
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Table  6.1  Binary  Interaction  Parameters  for 

the  NH3-HCI-H2O  System 

Binary  Pair  Bij 

Source 

Source 

0 

(e) 

0 

(e) 

h'*’-nhJ  0 

(e) 

0 

(e) 

0.1775 

(a) 

0.2945 

(a) 

0.20B 

(b) 

0.018  + 

(b.c) 

3.066?^ 

0.015 

3(aq) 

(b.c.d) 

0 

(b.c.d) 

nh^-nh!  0 

4  4 

(e) 

0 

(e) 

NH^-Cr  0.0522 

4 

(a) 

0.1918 

(a) 

NhJ-0H~  0.06 

(b) 

0.018  + 

(b.c) 

3.06B;j 

K-'‘”3(aq)  0.0117 

(b) 

-0.020 

(b) 

O 

1 

1 

o 

1 

o 

(e) 

0 

(e) 

o 

1 

X 

o 

1 

1 

(e) 

0 

(e) 

Cr-NH-_„,  0 

3(aq) 

N/A 

0 

N/A 

oh'-oh'  0 

(e) 

0 

(e) 

0»“-''»3(aq)  0-227 

(b.c) 

0 

(b.c) 

-  1.47E-3T 

4-  2.6E-6t2 

'*»3(aq)-N»3(aq)  '  0-026 

(b,c) 

0 

(b.c) 

+12.29/T 

(a)  Pitzer,  1979 

(b)  Maurer,  19B0 

(c)  Edwards  et  al . ,  1975 

(d)  Edwards  et  al . ,  197B 

(e)  Interactions  between  like-charged 

ions  are  neglected. 
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In  order  to  determine  the  accuracy  of  Edwards'  formulation  for  the 
species  under  consideration  In  this  study,  a  comparison  was  made  of  the 
theoretically  derived  and  experimentally  measured  activity  coefficients  for 
two  strong  electrolytes.  Figures  6.1  and  6.2  show  a  comparison  of  activity 
coefficients  calculated  using  Edwards  model  with  experimental  data  for 


hydrochloric  acid  at  25*  C  and  50*  C  respectively.  Figure  6.3  shows  the  same 
comparison  for  aqueous  solutions  of  at  25*  C.  In  all  cases  the 

agreement  Is  extremely  good  at  low  to  medium  concentrations,  but  significant 
deviations  are  present  at  high  concentrations.  This  may  be  due  to  the  Fact 
that  the  binary  Interaction  parameters  for  these  strong  electrolytes  are 
generally  fit  to  experimental  data  of  no  greater  than  6  molal  concentration. 
Predicted  results  exhibiting  similar  deviation  at  high  concentration  for  the 
activity  coefficient  of  hydrochloric  acid  have  been  obtained  by  Zemaltls 
(1986)  using  several  state-of-the-art  theoretical  models.  Including  Pitzers' 
original  formulation.  Clearly  the  art  of  modeling  strong  electrolytes  is  not 
developed  to  the  point  where  solutions  at  high  concentrations  can  be  predicted 
accurately. 

6.4  Activity  Coefficients  In  the  Gas  Phase 


The  gas  phase  was  assumed  to  be  non-ideal  due  to  the  highly  polar  nature 
of  the  molecules  present.  As  suggested  by  Prausnitz  et  al.  (1980)  the  virlal 
equation  was  chosen  as  the  mechanical  equation  of  state.  The  virial  equation 
obtained  from  a  power  series  expansion  of  pressure  can  be  written  as: 


Pv/RT  =  1  +  B  P/RT  +  C  (P/RT)‘  +  ... 
m  m 


(6.4) 


where  v  is  the  molar  volume  of  the  mixture,  and  B„  and  C„  are  the  second 

m  m 

and  third  virial  coefficients  for  the  gaseous  mixture,  respectively.  At  low 
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to  moderate  densities,  the  virial  equation  may  be  truncated  after  the  second 
term; 


Pv/RT  =  1  +  B  P/RT  (6.5) 

m 

According  to  Prausnitz  et  ai.  (1980,  p.  28)  the  virial  coefficient  for  a 
mixture  may  be  calculated  from  the  individual  virial  coefficients  of  the 
components  by  the  following  relation; 


B„  (T,  X 


m 


1’  2’ 


^ngas^ 


ngas 
=  I 

i=l 


ngas 

I 

j=l 


X.X.Bi. 


(T) 


(6.6) 


The  B^j  are  the  individual  pure  component  (i  equal  j)  and  cross  component  (i 
not  equal  j)  second  virial  coefficients.  The  activity  coefficients  may  be 
obtained  from; 


In 


ngas 


(5.7) 


The  method  of  Hayden  and  O'Connell  (1975)  was  used  to  calculate  the 

B.j.  The  necessary  input  parameters  for  each  component  are;  critical 

temperature  T,,  critical  pressure  P^,  dipole  moment  D*,,  and  mean  radius 
C  c  n 

of  gyration  Rp,  In  addition,  an  association  parameter  is  needed  for 
each  binary  pair.  These  parameters  are  given  for  a  large  number  of  species  by 
Prausnitz  et  al.  (1980)  pp.  145-178.  The  values  used  in  this  study  are 
presented  in  Tables  6.2  and  6.3.  Prausnitz  also  presents  the  numerous 
equations  (approximately  30)  necessary  to  calculate  the  B^  (pp.  130-133). 

They  will  not  be  repeated  here. 
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Table  6.2.  Molecular  Species  Data  for  Virial  Model 


Species 

Tc 

(K) 

Pc 

(bar) 

Ro 

(angstrom) 

Om 

(Debye) 

H^O 

647.37 

221 .20 

0.615 

1 .83 

NH3 

405.54 

112.80 

0.853 

1 .47 

HCl 

324.54 

82.60 

0.299 

1 .07 

Table  6.3. 

Binary  Association 

Parameter 

j  for  Virial  Model 

“20(g) 

““3(g) 

HCl,g, 

*^20(9) 

1.70 

0.0 

1  .38 

NH3(g) 

0.20 

0.0 

2.20 

HCl(g) 

0.0 

0.0 

0.0 
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In  order  to  check  the  accuracy  of  the  model,  the  virial  coefficient  of 
steam  was  calculated  using  Hayden  and  O'Connells'  method  (Figure  6.4).  The 
results  are  identical  with  those  presented  by  Prausnitz  et  al.  (1980),  p.  30. 

In  addition,  the  specific  volume  of  steam  was  calculated  from  the  virial 
equation  as  follows 

V  =  RT/P  +  (6.8) 

m 

These  results,  along  with  results  using  the  ideal  gas  equation,  are 
compared  against  experimental  results  in  Figure  6.5.  The  agreement  near  the 
boiling  temperature  is  noticeably  better  with  the  virial  equation. 


6.4  Binary  System  Equilibrium  Calculations 

Using  these  activity  coefficient  relationships,  the  computer  program  was 
tested  for  three  binary  systems  in  which  experimental  data  was  available. 
Figure  6.6  shows  predicted  and  experimental  solubility  data  for  the 
NH^Cl-H^O  system.  Agreement  is  good  at  lower  temperatures,  but  deviates 
significantly  at  higher  temperatures.  It  is  felt  that  this  deviation  at  high 
temperatures  may  be  due  to  the  variation  of  the  partial  molar  specific  heat  of 
NH^  and  Cl"  with  temperature.  Values  for  these  ions  were  only 
available  at  298.15  K.  Figures  6.7  and  6.8  show  predicted  and  experimental 
boiling  diagrams  for  the  system  and  for  the  HCl-HpO  system, 

respectively.  It  can  be  seen  that  the  agreement  for  the  weak  electrolyte, 
NH3(aq),  is  excellent,  while  the  agreement  for  the  strong  electrolyte, 

is  only  fair.  This  discrepancy  is  a  direct  result  of  the  inability 
of  the  activity  coefficient  model  for  HCl^^^^  to  accurately  predict  solution 
behavior  at  extremely  high  concentrations. 
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6.6  EQuilibrium  Calculations  for  HCI-NHj-H^O  System 

Basically  encouraged  by  the  test  cases,  the  computer  program  was  used  to 
predict  the  equilibrium  compositions  of  the  complete  HCl-NH^-H^O  ternary 
system.  The  version  of  the  computer  program  used  for  this  calculation  is 
listed  in  Appendix  A.  Results  for  0%,  10%,  20%,  and  30%  concentrations  of 
ammonia  in  water  are  shown  in  Figures  (6.9-16).  Mixture  density,  temperature, 
phase  and  species  mass  fractions  and  the  void  fraction  (volume  fraction  of  the 
gas  phase)  are  plotted  against  mixture  fraction,  f.  Calculations  presented  in 
this  form  (using  mixture  fraction  to  specify  the  relative  amount  of  reactants) 
can  be  used  directly  in  the  reaction  zone  modeling  technique  used  by  this  lab 
(Chan  et  a1.,  1987-1988).  Overall,  these  results  are  fairly  close  to  those 
obtained  by  previous  methods  for  this  system  (Chan  et  al.,  1987-1988).  The 
primary  difference  is  that  the  modeling  techniques  described  in  this  report 
were  able  to  predict  the  formation  of  a  solid  phase,  NH^Cl^^j. 
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MEHN  nCTIUITY  COEFFICIENT 


FIGURES.  1  riEPlN  PiCTIUITY  COEFFICIENT.  HCL  PT  25  DEB.  C 
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MEAN  ACTIUITY  COEFFICIENT 


FIGURE  6.2  riERN  PICTIUITY  COEFFICIENT,  HCL  AT  50  DEB.  C 


MEAN  ACTIUITY  COEFFICIENT 


FIGURE  6.2  MERN  RCTIUITV  COEFFICIENT,  NH4CL  PT  25  DEG.  C 


MOLALITY 
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(  CC/MOLE  ) 


FIGURE  6.4  SECOND  UIRIAL  COEFFICIENT  FOR  STEAM 
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TEMPERATURE  (C) 
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SPECIFIC  UOLUME  ( M**3/KB  ) 


FIGURE  6.5  SPECIFIC  UOLUtiE  OF  STEPM  AT  0.  1  MPA 


TEMPERATURE  (C) 
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TEMPERATURE  (C) 
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FIGURE  6.9  PURE  H20(  L)  +  HCLCG).  T I NF  =T0  =298K .  1  ftTM 
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TEMPERATURE  (K)  MASS  AND  UOID  FRACTION 


FIGURE  G.  10  PURE  H20(  L  )  >  HCL(  G  ).  TINF=T0=298K.  1  ATM 
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mss  AND  UOID  FRACTION  DENSITY  (  KG/M**3 ) 


FIGURE  6.11  10/  NH3C  AO  )  HCL(  G  ).  TINF*T0=29SK.  1  PTH 
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TEMPERATURE  (K)  MASS  AND  UOID  FRACTION 


■IGURE  6.12  IBZ  NII3(00)  +  IICL<6).  T  1NF=T0=298K,  1  OTM 
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FIGURE  B.14  20J!'.  NH3<  AO )  +  HCLCG).  T1NF=T0=298K .  1  ATU 
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MASS  AND  UOID  FRACTION  DENSITY  ( KB/Mf ¥3  ) 


FIGURE  6.15  30/  NH3(  AD  )  +  HCLC  G  ).  T I NF  =T0  =298K ,  1  ATM 


1000. 

100. 

10. 

1 . 

.  1 

1. 

.  1 

.  01 

.  001 


500. 
450. 
400. 
250. 
200. 

MIXTURE  FRACTION  -  F 


H 

m 

"D 

m 

i> 

H 

C 

73 

m 


7s 

'w' 


-  58  - 


FIGURE  G.  16  3Qx.  MIKK  DO  )  +  IICL(G).  TINF=T0=298K.  1  ATM 
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7.  APPLICATION  TO  THE  SOOIUM-CHLORINE  SYSTEM 


In  recent  years,  several  experimental  and  theoretical  studies  of  the 
structure  of  a  gaseous  chlorine  jet  discharging  into  a  stagnant  bath  of  liquid 
sodium  have  been  published  (Avery,  1974,  Chen  and  Faeth,  1982,  and  Cho  et  al., 
1967).  Present  work  by  this  lab  has  focused  on  extending  the  modeling 
techniques  of  Chen  and  Faeth  (1982)  to  include  the  effect  of  radiation  and 
enclosure  on  the  reaction  zone  (Chan  et  al.,  to  be  published).  A  necessary 
input  for  these  calculations  is  the  equilibrium  state  reached  after  mixing 
various  amounts  of  nozzle  and  bath  fluids.  In  this  system  the  bath  fluid 
(«)  consists  of  pure  liquid  sodium.  The  nozzle  fluid  (0)  is  pure 

chlorine  gas,  CI2. 

7.1  Gas  Phase 

Following  the  lead  of  Avery  (1974),  the  gas  phase  was  assumed  to  consist 
of  Cl^gj,  Cl^^gj.  Na2(g)’  "^^^^g)’  ^®2‘'^2(g)  ®  perfect  gas 

mixture.  Thus,  the  chemical  potentials  of  the  gaseous  species  were  described 
using  equation  3.25.  Ionization  of  the  gases  was  neglected,  as  preliminary 
calculations  revealed  extremely  low  concentrations  of  ionic  species,  even  at 
high  temperatures. 

7.2  Solid  Phases 

Both  the  metal,  Na^^^,  and  the  salt,  NaCl^^^,  were  considered  in 
single  species  solid  phases.  Consequently,  equation  3.33  was  used  to  describe 
the  chemical  potential  of  these  species.  Of  course,  the  solid  metal  does  not 
appear  in  the  system  unless  the  bath  temperature  (t“)  is  below  the  melting 
point  of  sodium  (371K).  The  bath  temperature  in  the  present  calculations  is 
above  the  melting  point  of  the  metal. 
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7.3  Liquid  Phase 

The  liquid  phase  consists  of  the  liquid  metal,  and  the  molten 

salt,  NaCl^^^j.  Bredig  (1964)  has  shovm  that  these  species  may  form  a  two 
phase,  immiscible  liquid  mixture.  The  lighter  liquid  is  sodium-rich, 
containing  a  small  amount  of  dissolved  salt.  The  heavier  liquid  is  salt-rich, 
containing  a  small  amount  of  dissolved  sodium.  Avery  (1974)  has  used  the  van 
Laar  model  (Lewis  and  Randall,  1961,  pp.  287-288)  to  accurately  predict  the 
mutual  solubility  of  these  two  species.  For  a  binary  solution  (denoting 
species  Na^^^j  and  NaCl^j^j  as  1  and  2,  respectively),  the  activity 
coefficients  may  be  expressed  as: 


In 


In  Y^  = 


(7.1) 

(7.2) 


where 


=  N^b^/(N^b^  +  N2b2) 
Z2  =  N2b2/(N^b^  +  N2b2) 


(7.3) 

(7.4) 


The  preceding  equations  contain  three  empirical  parameters  which  must  be 
fit  to  experimental  solubility  data  (b^,  b2,  and  A^2)-  ^e  can  reduce 

t 

this  to  two  parameters  by  defining  8  =  b2/b^  and  A^2  =  ^i^i2'  equations 
for  the  activity  coefficients  then  become: 


In  Y^  =  A^2  [(BN2/(N^  +  BN2)]  /RT 


(7.5) 


In  Y2  =  BA^2  [N.,/(N^  +  BN2)]‘/RT 


(7.6) 


Avery  (1974)  has  fit  A^2  ^nd  B  to  experimental  data, 
expressions  are 


The  resulting 
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B  =  0.929 


(7.7) 


A.,2  =  23.36  -  0.01325  T  (7.8) 

I 

where  units  of  Kcal/gmole  and  T  in  units  of  K. 

7.4  Thermodynamic  Data 

The  necessary  thermodynamic  data  for  this  system  was  obtained  from  data 
files  provided  with  the  CEC-72  computer  code  (Gordon  and  McBride,  1971). 

These  files  consist  of  the  polynomial  constants  for  curve-fits  of  the 
thermochemical  data  of  each  species.  The  original  source  of  the  data  for  most 
of  the  species  is  the  JANAF  Thermochemical  Tables  (1971).  Each  species  has  a 
set  of  14  coefficients.  Seven  coefficients  apply  to  the  low  temperature  range 
(300-1000  K).  and  seven  coefficients  apply  to  the  high  temperature  range 
(1000-5000  K) .  The  CEC-NMS  computer  program  was  modified  to  read  this  data 
file  and  calculate  the  specific  heat,  enthalpy,  entropy,  and  Gibbs  energy  of 
the  relavent  species.  For  a  listing  of  the  computer  code,  see  Appendix  B. 

7.5  Results 

The  equilibrium  results  for  T*  =  1130  K  nd  T®  =  298.15  K  are 
presented  in  Figures  7.1  and  7.2.  These  particular  conditions  were  chosen  to 
allow  direct  comparison  with  results  presented  by  Chen  and  Faeth  (1982).  The 
sodium-rich  liquid  is  denoted  by  L-1  and  the  salt-rich  liquid  is  denoted  by 
L-2.  There  is  essentially  no  difference  in  the  results  calculated  by  the  two 
different  methods. 
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HASS  AND  UOID  FRACTION  DENS I TV  ( KG/M**3  ) 


FIGURE  7.1  N«(L)-CL2(G)  SYSTEM.  T1NF=1130K.  T0=298K.  1  ATM 
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TEMPERATURE  (K)  MASS  AND  UOID  FRACTION 


FIGURE  7.2P  NA(  L  )-CL2(  G )  SYSTEM.  TINF  =  1130K.  T0=298K,  I  ATM 
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FIGURE  7.2B  NA(  L)-CL2< G )  SYSTEM.  TINF  =  113QK.  T0=298K.  1  HTM 
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8.  APPLICATION  TO  THE  LITHIUM-FLUORINE  SYSTEM 


A  system  similar  to  the  sodium-chlorine  reactant  pair  described  in  the 
previous  section  results  from  the  injection  of  gaseous  fluorine,  into 

a  stagnant  liquid  lithium  (Li^L)^  bath.  Once  again,  the  approach  used  is 
based  on  the  thermodynamic  work  of  Avery  (1974). 

8.1  System 

The  gaseous  phase  was  assumed  to  consist  of  *^2(g)’  ^\g)’  ^^2(g)’ 

*■^^(9)'  *’^2^2(g)’  '■^3^3(g)  ®  perfect  gas  mixture.  Ionization  was 

neglected  as  preliminary  calculations  showed  that  the  maximum  concentration  of 
ions  at  the  highest  temperatures  was  approximately  IX  by  mass. 

Li^^j  and  LiF^^^,  the  solid  metal  and  salt,  were  considered  as  single 
species  phases.  Immiscible  liquid  phases  consisting  of  dissolved  in 

Li^l^j  and  dissolved  in  LiF^j^^  were  modeled  as  before,  using  the 

I 

van  Laar  model.  The  coefficients  A^^  ®  were  provided  by  Avery  (1974) 

as 

8  =  0.725  (8.1) 

A^2  =  24.6  -  0.01072  T  (8.1) 

I 

where  A^j^  in  units  of  Kcal/gmole.  Note  that  is  designated  by 

1  and  sodium  case,  thermodynamic  data  was  read  in 

from  the  CEC-72  data  file. 

8.2  Results 

00  0 

The  results  of  the  equilibrium  calculations  for  T  =  1130  K  and  T 
=  298.15  K  at  1  atmosphere  are  shown  in  Figures  8.1  and  8.2.  No  calculations 
of  this  type  (based  on  mixture  fraction)  for  this  particular  system  have  been 
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presented  in  the  literature.  Consequently,  comparison  with  other  calculation 
methods  is  not  possible.  We  can  see,  however,  the  similarities  of  the 
lithium-fluorine  system  to  the  sodium-chlorine  system  described  in  Section  7. 
The  primary  differences  are  the  higher  boiling  point  of  the  liquid  lithium 
(1620  K  vs  1156  K  for  sodium),  and  the  higher  peak  equilibrium  temperature 
(4573  K  vs.  3352  K  for  sodium-chlorine  reaction). 

For  the  initial  conditions  presented  here,  the  liquid  phase  did  not  split 
into  two  immiscible  liquids.  For  these  calculations,  the  salt-rich  liquid 
does  not  form  until  the  temperature  of  the  mixture  is  well  above  the  consolute 
temperature  of  1603  K  (Bredig,  1964,  p.  373).  Note  that  the  consolute 
temperature  is  the  highest  temperature  that  the  two  liquids  coexist.  Had  the 
initial  bath  temperature  been  lower,  it  is  likely  that  the  solubility  limit  of 
the  salt  in  the  metal  would  have  been  reached  at  a  lower  temperature.  Then 
the  salt-rich  liquid  would  begin  to  form  a"-  a  temperature  below  the  consolute 
temperature.  In  this  case  an  immiscible  liquid  region  would  appear  as  in  the 
sodium  system  described  in  the  previous  chapter. 
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FIGURES. 1  LICL)-F2CG>  SYSTEM.  TINF=1120K.  T0=29SK.  1  ATM 
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TEMPERATURE  (K)  MASS  AND  UOID  FRACTION 


FIGURE  8.2  LI(L)-F2CG)  SYSTEM.  TINF  =  1130K.  TQ=298K,  1  ATM 
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9.  APPLICATION  TO  THE  LITHIUM-SULFUR  HEXAFLUORIDE  SYSTEM 


The  system  resulting  from  the  injection  of  gaseous  sulfur  hexafluoride 
(SFg(g))  into  liquid  lithium  (Li^^^)  is  discussed  in  this  section.  This 
particular  pair  of  reactants  is  under  consideration  as  a  thermal  energy  source 
for  underwater  applications  (Mattavi  et  al.,  1969;  van  der  Sluys,  1975; 
Biermann,  1975;  and  Groff  and  Faeth,  1978A).  The  thermodynamic  model  used  is 
based  on  the  work  of  Groff  and  Faeth  (1976,  and  1978B). 


9.1  System 

The  gas  phase  was  assumed  to  consist  of  the  species  Li^gj,  Li2^gj. 

^^*^(9)’  '■^2*^2{g)’  ^S^3(g)’  ^^6(g)'  *’^2^(g) 

gas  mixture.  These  preliminary  calculations  do  not  include  any  of  the  lower 
sulfur  fluorides  such  as  SF^,  etc.  It  is  likely  that  they  would  be  present 
in  significant  amounts,  especially  at  higher  temperatures.  The  presence  of 
ions  in  the  gas  phase  was  also  neglected.  Three  solid  species  were  considered 
in  separate  phases:  ^’2^(s)’ 

As  in  the  previous  two  sections,  the  liquid  metal,  and  molten 

salt,  LiF^j^j,  form  two  immiscible  phases.  In  the  sulfur  hexafluoride- 
lithium  system,  however,  there  is  an  additional  component,  LipS^l^j, 
dissolved  in  each  phase.  The  activity  coefficients  of  the  immiscible  liquid 
phases  are  predicted  using  the  van  Laar  model  extended  to  multicomponent 
systems  (Wohl,  1946).  The  activity  coefficient  of  the  kth  species  may  be 
given  as 


In  Y.,  = 


1  J 


RT  a  N.b.)' 
i  ’  ’ 


(9.1) 
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I  III 

where  A.,  and  b.  are  empirical  parameters.  Note  that  A..  =0,  and  A..  =  A... 
1j  1  n  ij 

The  equation  given  here  differs  from  that  presented  by  Groff  and  Faeth  (1976, 
p.  22,  and  1978,  p.  328).  We  believe  this  equation  to  be  the  correct  one. 

For  a  system  with  three  components,  the  number  of  empirical  parameters 
may  be  reduced  by  one  by  defining  new  parameters  A.j  and  B.  as  (Groff  and 
Faeth,  1976,  pp.  20-21): 


I 


(9.2) 


and 


B.|  =  b.|/b2;  82  =  b2/b2;  and  =  b^j/b^  =  B.|B2  (9.3) 

The  subscripts  1,  2,  and  3  refer  to  Li^L)’  '-’^(L)’  *-^2^’  'f'^^Pecti vely . 

Equation  9.1  may  then  be  written  as 


I  I  (A,,  -  0.5A,j)  (bj/b,) 


KT  cz  ",  Oi^l’k)!' 


(9.4) 


Groff  and  Faeth  (1978B)  have  estimated  the  empirical  parameters  from  experi¬ 
mental  data  as  a  linear  function  of  temperature  using  the  least  squares  tech¬ 
nique.  The  resulting  equations  are: 

A^2  =  121.53  -  0.05730  T  (kJ/mole)  (9.5) 

A^2  =  53.287  -  4.00E-4  T  (kJ/mole)  (9.6) 

A23  =  0  (9.7) 

and 

B^  =  2.5B46  -  8.132E-4  T  (9.8) 
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83  =  1.36 


(9.9) 


82  =  83/8^  (9.10) 

9.2  Thermodynamic  Data 

Thermodynamic  data  for  all  three  forms  of  Li2S  (lithium  sulfide)  is 
extremely  limited.  Groff  and  Faeth  (1976,  pp.  131-138)  estimated  the  required 
spectroscopic  data  for  determining  the  properties  of  Li2S^gj  using 
statistical  methods.  They  felt  the  accuracy  of  the  resulting  thermodynamic 
dat  for  Li2S^gj  to  be  within  +5%.  Their  equations  are  used  in  this  work 
to  calculate  the  enthalpy  and  Gibbs  energy  of  gaseous  lithium  sulfide.  Groff 
and  Faeth  estimated  thermodynamic  data  for  Li2S^^j  and  Li2^(L) 

Na2S,  Li20,  and  Na20  using  Kopp's  rule  (Lewis  and  Randall,  1961,  pp. 

57-58).  Their  estimated  values  were  used  to  calculate  the  enthalpy  and  Gibbs 
energy  for  solid  and  liquid  lithium  sulfide.  The  thermodynamic  data  for  all 
other  species  were  taken  from  the  data  file  provided  with  the  CEC-72  computer 
program  (Gordon  and  Mc8ride,  1971). 

9.3  Results 

CO  Q 

The  equilibrium  results  for  T  =  1130K  and  T  =  298K  are  presented 
in  Figures  9.1  and  9.2.  No  experimental  data  is  available  to  compare  these 
results  to,  nor  have  any  calculations  of  this  type  appeared  in  the 
literature.  As  in  the  lithium-fluorine  calculations  presented  in  the  previous 
section,  the  initial  conditions  are  such  that  salt-rich  and  metal-rich  liquid 
phases  do  not  appear  simultaneously  at  any  mixture  fraction. 
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FIGURE  9.1  LICL)-SF6(G)  SYSTEM.  TINF=1130K.  T0=298K. 
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10.  CONCLUSIONS 


The  computer  program  CEC-NMS  has  been  found  to  effectively  compute 
chemical  reaction  equilibria  for  a  variety  of  complex  systems.  Use  of  the 
direct  minimization  technique  allows  equilibrium  to  be  calculated  without 
regard  to  the  actual  reaction  equations  and  allows  easier  removal  of  trace 
species. 

For  the  majority  of  calculation  points,  convergence  was  achieved  in  under 
10  seconds  on  a  Uni  vac  1100  mainframe  computer.  While  not  unreasonable, 
future  work  might  include  modifications  to  the  program  to  reduce  this 
calculation  time.  This  would  become  more  important  if  the  program  were  to  be 
used  on  a  microcomputer  or  if  extremely  complex  systems  were  considered. 

While  an  attempt  was  made  to  make  the  program  user-friendly,  improvements 
in  this  respect  are  always  possible.  This  might  include  altering  the 
structure  of  the  program  to  improve  readability  and  developing  more  detailed 
documentation.  Another  Improvement  to  the  program  would  be  adding  a  routine 
for  automatically  checking  and  modifying  the  species  coefficient  matrix  if  it 
becomes  singular  during  the  course  of  computation.  Finally,  while  the  program 
has  converged  with  rather  arbitrary  initial  guesses  for  all  attempted  systems, 
there  are  algorithms  available  for  automatically  providing  initial  estimates 
based  on  the  Gibbs  free  energy  of  the  species  (Smith  and  Missen,  1982,  pp. 
201-204,  and  van  Zeggeren  and  Storey,  1970,  pp.  128-132).  The  inclusion  of 
such  an  algorithm  may  be  useful  for  certain  systems. 
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APPENDIX  A-0  PROGRAM  LISTINGS 


Legible  copy  can  not  be  reproduced. 
Interested  readers  should  contact  the  author  (S. 


H.  Chan). 
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